Existence and stability of multiple solutions to the gap equation 
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We argue by way of examples that, as a nonlinear integral equation, the gap equation can and 
does possess many physically distinct solutions for the dressed-quark propagator. The examples 
are drawn from a class that is successful in describing a broad range of hadron physics observables. 
We apply the homotopy continuation method to each of our four exemplars and thereby find all 
solutions that exist within the interesting domains of light current-quark masses and interaction 
strengths; and simultaneously provide an explanation of the nature and number of the solutions, 
many of which may be associated with dynamical chiral symmetry breaking. Introducing a stability 
criterion based on the scalar and pseudoscalar susceptibilities we demonstrate, however, that for any 
nonzero current-quark mass only the regular Nambu solution of the gap equation is stable against 
perturbations. This guarantees that the existence of multiple solutions to the gap equation cannot 
complicate the description of phenomena in hadron physics. 

PACS numbers: 12.38.Aw, 12.38.Lg, 11.30.Rd, 11. 10. St 



I. INTRODUCTION 

Dynamical chiral symmetry breaking (DCSB) is a par- 
ticularly striking feature of the Standard Model, play- 
ing an important role in formation of the visible mat- 
ter in the Universe It is apparent in the existence 
of a strongly momentum-dependent chiral-limit dressed- 
quark mass function, M(p 2 ), which is obtained in solu- 
tions of models for QCD's gap equation thatprovide a 
realistic description of hadron properties [2|-|j], and in 
a sharp increase in M(p 2 ) at p 2 < 2 GeV 2 when the 
current-quark mass is nonzero but light. The latter is 
seen in both Dyson- Schwinger equation (DSE) studies 0- 
uf and numerical simulations of lattice-regularised QCD 

FlOl] . This behaviour of the mass function must be part 
of any treatment of continuum strong QCD that aims to 
be considered realistic. 

The gap equation is a nonlinear integral equation. The 
nonlinearity gives it the power to express nonperturba- 
tive phenomena; and also leads to the curiosity that 
the solution is not unique. Mathematically, this should 
have been anticipated. However, perhaps surprisingly, 
it has only been established relatively recently that on 
a bounded, measurable domain of non-negative current- 
quark mass, realistic models of the QCD's gap equa- 
tion simultaneously admit more than one nonequivalent 
DCSB solution and also distinct solutions that may un- 
ambiguously be connected with the realization of chiral 
symmetry in the Wigner mode [ll| - [l3j . This feature can 
potentially create problems; e.g., if the additional solu- 
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tions are physically realisable and therefore have a mea- 
surable impact on observables. Thus, amongst the ques- 
tions that must be answered is that of which solution 
or solutions of the gap equation should be employed in 
defining the kernel of the Bethe-Salpeter two-body prob- 
lem. Naturally, in addressing such questions, the first 
thing to ensure is that one has found all solutions to the 
gap equation. 

Herein we describe solutions of the gap equation ob- 
tained using two different interaction models flil . flFjj and 
two vertex Ansdtze. We focus on the interesting domain 
of light current-quark masses and a large domain of in- 
teraction strengths. Notably, we employ a numerical 
method [HI], novel in the study of DSEs, which deliv- 
ers all solutions of the gap equation. This enables us to 
chart the complete solution set domains for each of the 
four kernels we consider. We also discuss the important 
issue of stability for each of the various solutions and 
thereby address the question raised above. 

Section[H]providcs a brief review of the gap and Bcthc- 
Salpeter equations. It also explains that all model- 
independent content of the so-called "Mexican hat" po- 
tential is encoded in the behaviour of the scalar and 
pseudoscalar susceptibilities, which therefore provide a 
general tool for judging the stability of gap equation so- 
lutions, insofar as any given solution might represent a 
valid starting point in the computation of hadron prop- 
erties. Section [IIII introduces the kernel Ansdtze that we 
employ. They are not new but were chosen because they 
are capable of describing and unifying a wide range of 
meson and baryon properties. Section IIVI is an exten- 
sive presentation and discussion of the solutions to the 
gap equation, which explains: their classification; their 
evolution in response to changes in two control parame- 
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ters (current-quark mass and interaction strength); and 
the process of identifying stable solutions. Section [V] is 
a recapitulation with some additional comments. We ex- 
plain our numerical method for solving the gap equation 
in App.[£] 



II. GAP EQUATION 

Since DCSB is a phenomenon emerging from the strong 
physics of drcssed-quarks, it is often studied via QCD's 
gap equation!]] 



S-\p) =Z 2 (i 1 -p + m bm ) 



+Zi Jj*D^(p-q)— llJt S{q)—T v {q,p), (1) 

where: D^„ is the gluon propagator; T v , the quark-gluon 
vertex; J d , a symbol representing a Poincare invariant 
regularisation of the four-dimensional integral, with A 
the regularisation mass-scale; m bm (A), the current-quark 
bare mass; and Zi^{( 2 , A 2 ), respectively, the vertex and 
quark wave-function renormalisation constants, with £ 
the renormalisation point. We employ the renormalisa- 
tion procedures of Ref. [171 ] and the same renormalisation 
point, C = 19 GcV. 

The gap equation's solution is the drcssed-quark prop- 
agator, which is commonly written in one of three equiv- 
alent forms: 

S(p) = -i 7 • pa v (p 2 , C 2 ) + a s ( P 2 , C 2 ) , (2a) 
= l/{i 1 -pA(p 2 ,( 2 ) + B(p 2 ,( 2 )}, (2b) 
= Z(p 2 ,( 2 )/lH-p + M(p 2 )}. (2c) 

The mass function, M(p 2 ), is independent of the renor- 
malisation point; and the renormaliscd current-quark 
mass, 



m C = Z m (C, A)m bm (A) = Z^Z 2 m hn \ 



(3) 



where Z± is the renormalisation constant associated with 
the Lagrangian's mass-term. The renormalisation-group 
invariant current-quark mass may be inferred via 
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(4) 



where j m = 12/(33 — 2Nf) with Nf the number of active 
quark flavours. The chiral limit is 



rh = . 



(5) 



1 Wc use a Euclidean metric: {7^,7^} = 2<5 M „; 7 J = 7 M ; 75 



74717273, tr[757 M 7„7 p 7 (T ] = -4e MJ/pCT ; a^ v 
a i b i> and p v- timclikc =>- P 2 < 0. 



Chiral-limit QCD possesses a SU L (N f ) ® SU R (N f ) 
symmetry and thus separates into two non- 
communicating theories: one for left-handed quarks 
and another for right-handed quarks. This can be 
seen to entail that the regular parts of the scalar and 
pseudoscalar vacuum susceptibilities must be identical 
(la ]. In fact, this is the content of the so-called "Mexican 
hat" potential, which is commonly used in building 
models for QCD. 

The symmetry requires that the gap equation is in- 
variant under a change in the sign of B(p 2 ) in Eq. (|2bl) : 
i.e., if Bo is a solution, then so is (—Bo). In the con- 
text of simp le g ap equation truncations, this has long 
been known [ljj, [20] ■ On the other hand, as we now de- 
scribe, it limits what may be called realistic Ansatze for 
the dressed-quark-gluon vertex. The dressed vertex sep- 
arates into a sum of two pieces: one in which every term 
is even in the number of Dirac matrices, r^ _ovcn ; and 
another in which every term is odd, r^ _odd . The vertex 
satisfies its own DSE, the kernel of which features the 
dressed-quark propagator. The result described above 



entails that T 



D- 



obtained as a solution of this DSE 



must change sign under B — > —B but is otherwise un- 
changed, whereas r^~ odd is invariant under 
This is a significant result because the scalar functions 
that accompany the even and odd tensor structures in the 
vertex cannot in general be expected to share the func- 
tional form of cither A(p 2 ) or B(p 2 ), or simple combina- 
tions thereof. Notwithstanding this, most vertex models 
are constructed with a simple functional dependence on 
A(p 2 ) and B(p 2 ), such that they do exhibit the correct 
properties under B —B [2l| - (26| . 

In the chiral limit, therefore, at least two possibilities 
are realisable. Namely, if the support of the integrand in 
the equation for B(p 2 ) is too small, then the gap equa- 
tion possesses solely the B = Bw = solution, whereas, 
if the support exceeds some critical value, it has three 



solutions; viz., B\y, B 



B N 



Bo, B — B 



N 



B . In 



studies that convert the gap equation into a second-order 
nonlinear differential equation for B(p 2 ) [2?], a step 
which is quantitatively accurate for p 2 > 2GeV 2 , it is 
natural to characterise the latter two solutions as regu- 
lar, whereas the first solution can be connected with an 
irregular solution. It is more common, however, to de- 
nominate the B = solution as the Wigner-Weyl mode 
and the latter two as Nambu-Goldstone type solutions, 
since they signal the dynamical breaking of chiral sym- 
metry. 

It is notable that in the chiral limit the Nambu so- 
lutions are energetically favoured in concrete computa- 
tions that produce both the Wigner- and Nambu-type 
solutions [ioi • In the context of Ref. [3(| , the continuum 
of Nambu solutions would be described as equivalent de- 
generate vacua. If one introduces a small current-quark 
mass, then solutions smoothly connected to B^, B~^ and 
Bw persist pdl - [l3l ]. Plainly, therefore, mere existence as 
a solution of the gap equation does not guarantee that 
solution's stability. 
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FIG. 1. Classical potential often imagined in connection with 
dynamical chiral symmetry breaking. The point "I" is the 
global minimum, characterised by the conditions in Eqs. iJBJ; 
"S" is a saddle-point, m 2 a > but mi < 0; and U is an 
unstable local maximum m 2 < and m„ < 0. 



In investigating stability of solutions to the gap equa- 
tion it has proved useful to employ the chiral suscepti- 
bility, which is defined as usual via the scalar vacuum 
polarisation. A solution is energetically unstable in re- 
sponse to fluctuations of some source if the associated 
chiral susceptibility is negative [29|, [3ll - l33l |. The infor- 
mation is incomplete, however, since if the susceptibility 
is positive semi-definite, then the solution may be sta- 
ble, meta-stable, or a saddle-point. The last possibility 
is real here because the scalar and pseudoscalar vacuum 
polarisations are distinguishable when chiral symmetry 
is dynamically broken [18| , and the pseudoscalar suscep- 
tibility can be negative whilst the scalar susceptibility is 
positive semi-definite. 

These polarisations are obtained via second-order func- 
tional derivatives of the theory's generating functional for 
connected one-particle-irreducible Schwinger functions; 
viz., 5 2 Tipi/SJ(x)5 J(y), where J — S,P are respectively 
scalar and pseudoscalar sources. Their importance and 
relevance herein are evident once one appreciates that a 
simultaneous consideration of the scalar and pseudoscalar 
vacuum polarisations expresses, amongst other things, 
all model-independent physical content of the so-called 
"Mexican hat" potential [33| • Using this connection here 
for illustrative simplicity, suppose that potential is rep- 
resented as U[s,p]. In this case an extremum v = {s,p}, 
is stable if, and only if, 



1 92 m i 

2 d* U[8 > p] 

1 92 m i 

2 W U[SM 



> 0, 



s.p 



> 0. 



(6a) 
(6b) 



If just one of m 2 , m 2 is negative, then v is a saddle-point; 
whereas if both are negative, then v is a local maximum. 
These points are depicted in Fig. [I] 

Translating back to the general case, Eqs. ((SJ) corre- 



spond to the statement that a solution configuration is 
truly stable if, and only if, both susceptibilities are pos- 
itive at zero total momentum. One can extend this by 
noting that each susceptibility is the inverse of a fully- 
dressed propagator for composite correlations in the rel- 
evant channel. It follows that Eqs. ((6]) translate into the 
statement that stability is guaranteed if, and only if, the 
lowest mass excitation in each channel has positive mass- 
squared. This is the criterion we use subsequently It can 
be implemented simply by solving the inhomogeneous 
scalar and pseudoscalar Bethe-Salpeter equations [3^ |: 

,A 

Tj(q; K) = Z 4 Mj - / g 2 D^(q - t) 

Jdi 

x^s(£ + )rj(e;K)s(i-)^r u (i_, q _) 

+ £ g 2 D^(q - e)^ S (£ + )^-AMq, I- K) (7) 

wherein the dressed-quark propagator is that which char- 
acterises the gap equation solution whose stability is in 
question and, furthermore: M{g P } = {1,75}; £± = 
£±K/2, without loss of generality in our Poincare covari- 
ant approach, where K is the total momentum entering 
the vertex; Aj u is a Bethe-Salpeter kernel, which is fully 
determined by the kernel of the gap equation; and the 
Bcthc-Salpctcr amplitudes have the general form 

T s (q; K) = [E s (q; K) + ^ 7 • K F s (q; K) (8a) 
+ 17 ■ qG s (q;K) + Hs(k;P)}, 

T P (q; K) = l5 [iE P (q; K) + 7 • K F P (q; K) (8b) 
+ 7 • q G P (q; K) + a^q^K^Hp (q; K)] , 

with Ej, Fj, Gj, Hj being scalar functions. We locate 
the lowest-mass excitation using the method of Ref. (3|| , 
which simplifies computations by permitting one to em- 
ploy solely spacclike momenta. 



III. MODEL KERNELS 

The gap equation's kernel is specified by the form 
used to express the contraction Zig 2 Dfj, v (k — q)T„(q,p) 
in Eq. ([T]). Herein we compare four kernels, which may 
all be introduced by first writing (k = p — q) 

z l9 2 D^(k)r u ( q ,p) = k 2 g(k 2 )D^ c (k)r^(q, P ) 

= [k 2 Gm(k 2 ) +4™ pQCD (fc 2 )] D^(k)T^q 7 p), (9) 

wherein: D ! ™ c (k) is the free-gauge-boson propagator^ 
<ipQCD(fc 2 ) is a bounded, monotonically-decreasing reg- 



2 We use Landau gauge, a choice made for many reasons [2(1 139 . 
[37J, for example, it is: a fixed point of the renormalisation group; 
that gauge for which sensitivity to model-dependent differences 
between Ansdtze for the fermion-gauge-boson vertex are least 
noticeable; and a covariant gauge, which is readily implemented 
in simulations of lattice regularised QCD (see, e.g., Refs. [8l.l9ll38|- 
|43|| , and citations therein and thereto). 
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ular continuation of the pcrturbativc-QCD running cou- 
pling to all values of spacelike- fc 2 ; Gm(k 2 ) is an Ansatz 
for the interaction at infrared momenta: Qm(k 2 ) <C 
%>QCD(fc 2 ) Vfc 2 > 2GeV 2 ; and T^{q,p) is an Ansatz for 
that part of the dresscd-quark-gluon vertex which cannot 
be absorbed into Q{k 2 ). 
In all instances we use 1171 



47ra pQCD (fc 2 ) 



87r 2 7 m F(s) 



ln[r + (1 + s/A 2 ) 2 } 



(10) 



where: lm = 12/(33-2jV/), N f = 4, A QCD = 0.234 GeV; 
t = c 2 - 1; and J'(s) = {1 - exp(— s/[4m 2 ])}/s, m t = 
0.5 GeV. For the infrared, we compare two forms; viz., 
00 



Att 2 



-Dse 



(11a) 
(lib) 



These arc actually one-parameter models because in 
both cases there is a domain of u> throughout which, in 
rainbow-ladder truncation - see below, computed prop- 
erties of ground state vector and flavour-nonsinglet pseu- 
doscalar mesons 0, 0, , and nucleon and A proper- 
ties 0, 0] are almost unchanged along the trajectory 
Dlu = constant. That domain is us G [0.3, 0.5] GeV for the 
interaction in Eq. (|lla[ ), whereupon Dlj = (0.72 GeV) 3 
provides a good description of the observables identified 
0|; whilst for Eq. (|TTE|) the domain is a; € [0.4, 0.6] GeV, 
with Doj = (0.8 GeV) 3 providing the best achievable 
phenomenological results 0]. We will use the mid- 
point of each domain for computations throughout; i.e., 
w = 0.4 GeV in Eq. (fTTaj) and u = 0.5 GeV in Eq. ([Tib]) . 
Wc employ two simple models for the vertex: 



^lBC 



{q,p) = iv 



A(q 2 ) + A(p 2 ) 



(12a) 
(12b) 



The first implements a rainbow-ladder truncation of the 
gap and Bethe-Salpeter equations, which is the leading 
order in a widely-used symmetry-preserving DSE trun- 
cation scheme 0, 0] . The second model is a truncation 
of the Ball-Chiu Ansatz 211. It is far from the most 



general form [23|, |26[ but, in circumstances we expose, 
it produces some qualitative changes in Eq. ([1]) and thus 
serves to highlight the impact of a dressed vertex on the 
number and nature of solutions to the gap equation. 

In general one can construct the Bethe-Salpeter ker- 
nel, Aj;,(<7, £; K), associated with any T v (q,p) using the 
formulae in Ref. 0]. Herein, however, Aj v (q, £; K) is 
omitted for reasons we now explain. This term is iden- 
tically zero in rainbow-ladder truncation J34|; i.e., with 
Eq. p2ap . Hence the omission need only be discussed 
in connection with Eq. (|12b[) . Firstly, Eq. (|12b[) has the 
same Dirac structure as Eq. (|12a[) and hence the associ- 
ated Aj v {q, £; K ) cannot realistically have a large effect 
on masses obtained via the Bethe-Salpeter equation since 



it does vanish identically using Eq. (|12a[) . More generally, 
we use the Bethe-Salpeter equation primarily in order to 
gauge stability of solutions to the gap equation. A so- 
lution is stable if, and only if, both the scalar and pseu- 
doscalar mass-squared values are positive semi-definite 
when computed using that solution. In the scalar chan- 
nel the omission of As„(g,£; K) suppresses repulsion and 
hence produces a lower bound on the absolute value of 
the mass-squared 0. In the pseudoscalar channel, on 
the other hand, the diagrams represented by Ap„(q, £; K) 
largely cancel amongst themselves in the neighbourhood 
of the chiral limit, so this term has a negligible impact on 
the mass-squared within this domain [5ll . |52| . Hence the 
omission of Aj v (q, £; K) cannot materially affect a study 
of stability. 

It is worth remarking that, irrespective of the remarks 
just made, all kernels constructed using Eqs. (JTTJ) , ([12")) 
preserve the one-loop renormalisation-group behavior of 
QCD in the gap and Bethe-Salpeter equations. In the in- 
frared, on the other hand, there are differences between 
Eq. pTaj) and (|llbj) . They are detailed in Ref. 0; and 
notable amongst them is the fact that interactions con- 
structed from Eq. (|llbl) possess an infrared momentum- 
dependence that is consonant with modern DSE- and 
latticc-QCD results, whereas those produced by Eq. (|lla|) 
violate this constraint. Whilst this does not appear 
to impair the utility of Eq. (| 1 1 a.|) in connection with 
properties of ground state vector and flavour-nonsinglet 
pseudoscalar mesons, and the nucleon and A, it does 
markedly affect its predictions for quantities more sensi- 
tive to the infrared behaviour of the interaction, such as 
the properties of excited states 0, 0] and, as we shall 
see, the location of phase boundaries. 



IV. RESULTS AND DISCUSSION 

A. Solutions of the Quark DSE 

In solving the gap equation we have two control pa- 
rameters upon which the existence and number of solu- 
tions will depend: current-quark mass, m; and interac- 
tion strength, which will hereafter be characterised by 
the dimensionless numbe 101 = D/oj 2 . That DCSB is a 
possibility for to = and X > I c , where I c is some crit- 
ical value, guarantees that the gap equation does admit 
more than one solution: I c is a rh = bifurcation point 
0, . With the existence of furcation points assured, 
one must anticipate that the straightforward iteration 
procedure used commonly to solve the gap equation will 
be inadequate to the task of locating all solutions and 



3 For reference, in rainbow-ladder truncation the domain of rea- 
sonable interaction strengths; i.e., strengths with which one can 
hope to describe hadron phenomena, is 3 < X < 13 for Eq. Illlal 
and 2 < X < 8 for Eq. HTbll . 
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FIG. 2. Solutions of the gap equation obtained using mr = 
5MeV with X — 5.8 for the rainbow-ladder vertex, Eq. (|12a|l . 
or I = 3.1 for the 1BC vertex, Eq. ()12b|l . Upper panels: 
interaction of Eqs. (JTUJ, (|lla|l . {T5J); lower panels: Eqs. (fTUl) . 
()llb|) . (|12[) . All panels: solid curve, positive Nambu solution 
- N + ; dashed, Wigner solution; and dotted, negative Nambu 
solution - N~ . The insets highlight the infrared behaviour of 
the Wigner and A~ solutions, in particular their single zero. 



tracking their evolution as {rh,X} are varied. In con- 
trast, the homotopy continuation method, summarised 
in App.[5] is well suited to this challenge. 
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1. Influence of interaction strength 

To illustrate the point and establish a context we 
solved the gap equation with the four interaction ker- 
nels described above. The resulting mass functions are 
depicted in Fig. [2 With the listed parameters, the gap 
equation possesses three distinct solutions, as elucidated 
in Refs. [Ill4l3l ]. The figure displays one novelty, however: 
viz., both interactions support three nontrivial solutions 
with a dressed-vertex. This dressing, albeit apparently 
simple, does qualitatively change the gap equation, as we 
will explain below. 

In order to determine the solution set of the gap equa- 
tion, which, recall, is a pair of coupled, nonlinear inte- 
gral equations for two functions, we solved Eq. (I) on a 
large domain of {rh,I} £ M 2 . The values and parameter- 
dependence of the computed quantities A(0) and B(0) 
are useful in characterising the solutions. Some of this 
information is portrayed in Fig.[3J It is immediately ap- 
parent that, in the chiral limit, three critical points exist 
within the domain displayed. 



FIG. 3. X = D/oj 2 -dependence of A(0) and B(0) obtained 
using the models specified by Eqs. (|10[) . (|12a[) : upper 

grouping - Eq. (|llal) : and lower grouping - Eq. ()llb|) . The 
upper two panels in each grouping were computed in the chiral 
limit, whereas the lower two panels were obtained with mr = 
5 MeV or m" = 0.5 MeV, respectively. All panels: solid curve 
- N + solution; long-dashed - N~ ; short-dashed - regular 
Wigner solution; and dotted - a second Wigner-type solution. 
Naturally, in the chiral limit A N + = A N - ; and deviations 
from this identity are almost imperceptible for rh nonzero 
but small. 



The first is a trifurcation point. For X < If, 
the magnitude of which depends on the form chosen 
from Eq. pip , only the long known chiral-symmetry- 
preserving (Wigner) solution is present, which we here- 
after denote W or W\. At I = If, however, two new 
solutions appear. These are the normal DCSB (Nambu) 
solutions, described above, which we henceforth denote 
N + or and N~ or N± , respectively. 
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FIG. 4. Momentum dependence of the W2 solution obtained 
using Eqs. (TO}, ifTTa)) . (fl2a)) with 1 = 5.8 and m c = 5MeV. 
N.B. M(p 2 ) = for 171 — 0. 
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FIG. 5. I = D/cj 2 -dependence of A( 0) an d B(0) obtained 
using the model specified by Eqs. (|fOI) . ()IIb|l . (|I2a|l and with 
vnr — 3MeV. Solid curve - N + solution; long-dashed - N~ ; 
short-dashed - regular Wigner solution; and dotted - a second 
Wigner-type solution. 



The second critical point is associated with the appear- 
ance of a novel solution first observed in Ref. At 
X > 25 1 again interaction-dependent, a second Wigner- 
like solution, W2, appears: whilst Bw 2 {p 2 ) = 0, Aw 2 {p 2 ) 
is nontrivial and Aw 2 (p 2 ) ^ Ay/^p 2 ). The momentum- 
dependence of the W2 solution is depicted in Fig. [4] for 
small but nonzero current-quark mass. Plainly, the mass 
function has two zeros. 

A third critical point I = 1% locates an interaction 
strength at which the W\ and W2 solutions merge and 
beyond which they disappear. In the context of App.lAl 
it is a turning point. This is a novel result; for whilst the 
equation for B(p 2 ) always admits the B = solution, the 
existence of If indicates that if the coupling strength is 
strong enough, then the equation for A does not possess 
a solution. We have thus exposed two chiral-limit exam- 
ples of gap equations that only support a nonperturba- 
tive chiral symmetry preserving solution on a bounded 
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FIG. 6. I = D/cj 2 -dependence of ,4(0) and M(0) ob- 
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The upper two panels in each grouping were computed in 
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long-dashed - N~; short-dashed - dot-dashed - 7V^~; 

and dotted - Wigner solution. Naturally, in the chiral limit 
^4^+ = A N -, i — 1,2; and deviations from these identities 
are almost imperceptible for rh nonzero but small. 



domain of interaction strength. Hence, one cannot in fu- 
ture assume that a gap equation will always admit a fully 
self-consistent Wigner solution at strong coupling. 

The picture changes somewhat at nonzero current- 
quark mass. Of particular note: whilst the solution 
is always present, the W\ and TV-f solutions exist only on 
a domain I > If m > If: I = I° m is a turning point. 
Figure [2] shows that in this case Mw{p 2 ) is nonzero and 
both Mw(p 2 ), M N -(p 2 ) possess a zero. Also striking is 
the sensitivity of the Wi, W2 solutions to the infrared 
behaviour of the interaction, which is evident via com- 
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parison of Figs.0] and [5] In the neighbourhood of the 
chiral limit, both Wi and W% are absent for X > If, irre- 
spective of the interaction. Note, however, that one must 
be very close to rh = for this to be true when the inter- 
action is constructed using Eq. (|llb[) . In this case there 
is a current-quark mass, mj? , above which both W\ and 
W2 survive and evolve smoothly on X > X^ m > I| (see 
Fig.[5j). This is actually also true when the interaction is 
constructed using Eq. (|lla[) but m^ T /m^ c > 10. 

The preceding few paragraphs described properties of 
solutions obtained with gap equation interaction kernels 
built in the rainbow truncation; i.e., using Eq. (|12a| . 
Results obtained with the modestly dressed vertex in 
Eq. (|12bl) arc displayed in Fig. [6] A comparison of Fig. [6] 
with Figs. G2[5] reveals a significant difference; viz., using 
Eq. (|12b[) there is only ever one Wigner-type solution. A 
little algebra explains why: using Eq. (|1 2b[) the equation 
for A{p 2 ) derived from Eq. ((1) is actually a linear equa- 
tion for ay(p 2 ) when B = 0; and linear equations have 
at most one solution. (N.B. This is also true when the 
"2BC" vertex is used; i.e., the vertex obtained from that 
in Ref. [2l| by dropping only the Dirac-scalar term.) 

On the other hand, with increasing interaction 
strength, the number of Nambu-like solutions also grows. 
The solutions we've labelled as N^ 1 each possess a single 
zero in the chiral limit, irrespective of the choice made in 
Eq. ; and the 7V^~ solution has two zeros when > 
as a result of being required to equal this positive mass at 
the renormalisation point (see Fig. [7]). The momentum- 
dependence of the new Nambu solutions becomes quite 
complicated as the interaction strength reaches large val- 
ues. Notwithstanding this, there is always a value of the 
interaction strength above which these solutions exhibit 
the hallmarks of the normal Nambu solutions; i.e., in the 



FIG. 8. Current-quark-mass-dependence of A(0) and B(0) 
obtained using the model specified by Eqs. {TU|, fTTa}, (|12al) : 
upper pair - I = 3.5; middle pair - I = 5.5; and bottom pair 
- 1 = 7.8. All panels: solid curve - N + solution; short-dashed 
- Wi; long-dashed - N~; and dotted - W2. 



chiral limit they are nonzero mirror image pairs, and for 
small nonzero current-quark mass the members of the 
pair have commensurate magnitudes. 

We define Nambu-like to mean solutions with a 
nonzero mass function in the chiral limit. However, our 
nomenclature is not without ambiguity. For example, on 
1 < X < 6 the N2 solution has properties characteristic 
of a Wigner solution: the mass function is zero and it tri- 
furcates from the regular Nambu solutions at the lower 
boundary of this domain, evolving thereafter within the 
domain as a chirally symmetric solution. At the upper 
end, however, it trifurcates instead as the partner to the 
N% solution and maintains that DCSB trajectory. It is 
finally for this reason that we label it a Nambu-like solu- 
tion. This pattern might repeat again with increasing X, 
so that the solution we have labelled as Wigner-like is, in 
fact, the N3 solution, which will, in turn, trifurcate to 
form the DCSB partner of the N$ solution, leaving either 
a true Wigner solution or a N± solution, if the pattern 
is interminable. Whilst this is of academic interest it is 
not physically relevant since the values of X involved far 
exceed the upper bound on values which arc capable of 
producing an efficacious hadron physics phenomenology 
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FIG. 9. Current-quark-mass-dependence of yl(0), B(0) ob- 
tained using the model specified by Eqs. (| 10f) . (|llbl) . (|12al) : 
upper pair - X = 3.5; middle pair - X = 6.4; and bottom pair 
- X = 7.8. All panels: solid curve - N + solution; short-dashed 
- Wi; long-dashed - iV~; and dotted - W2. 



2. Influence of current-quark mass 



It will already be plain from Sec. II V All that the nature 
and number of solutions to the gap equation also depend 
on the current-quark mass. This is emphasised by Figs. [51 
[51 which show that the simultaneous existence of distinct 
solutions depends sensitively on the location in R 2 of the 
control parameters {to, I}. The Wigner solutions are 
again a good example. There are points {rh,X} G M 2 
at which A\y 2 (0) = and only with sufficiently large 
interaction strength is there a clear relationship between 
the Wi and W2 solutions. At such strengths, however, 
the Wigner solutions do not exist in a sizeable connected 
domain containing the chiral limit. (These points were 
mentioned earlier, in connection with Figs.[3j[5]) 

The influence on the solutions of dressing the quark- 
gluon vertex is illustrated in Fig. llQI In important re- 
spects, the picture is simpler in this case. As usual, the 
regular Nambu solution is distinct but all other solutions 
can be considered to appear in pairs, something we noted 
earlier in connection with Fig.|7J] One simple observa- 
tion is important and supported by the figure; viz., at 
fixed interaction strength the number of solution pairs 
decrements uniformly as the current-quark mass passes 
discrete critical values until only the regular Nambu so- 
lution exists. 



< 



L MT+1BC 




- 

_ J"-) =7.81 ; 

1 1 I 1 1 I 1 1 I 


— i 
✓ 

✓ 

: w n 2 

- 1 1 1 1 1 1 1 1 1 


— N , 

:/ ) 
-/ „ " 




D/w 2 - 
= 15.6 - 


£ - = 

-\ 

— 1 1 1 1 1 1 1 1 1 





1 



-1 o 
2 Q 

CD 
< 







0.00 0.03 0.06 0.09 0.03 0.06 0.09 

m [GeV] 



- -2 



< 





QC+1BC 


















/ 










/ 

/ 






D/w 2 




_ / 
_/ 


— n; 




=4.0 




- 1 1 1 1 


— N 2 

■ 1 1 1 I 1 1 


■ 1 1 1 1 


1 1 1 1 1 


1 : 




^1 










- -n; 








* S. 






D/w 2 




" J ^> 


, — w 

\ 




=16.0 














- 1 1 I I I 1 1 I I 1 I 1 


■ 1 1 1 1 


■ 1 1 1 1 





1 



A -1 



CD 
2 ^ 



0.00 0.05 0.10 0.05 0.10 

m [GeV] 

FIG. 10. Current-quark-mass-dependence of ^4(0), M(0) ob- 
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curve - N± solution; long-dashed - iV-f ; short-dashed - 7V^" ; 
dot-dashed - jV 2 ~; and dotted - Wigner solution. 



B. Solution Set Domains 



The results described in Secs. lIV A llllV A 21 ind icate 
that while choosing between Eqs. (|lla[) and (| 1 lh>|) pro- 
duces quantitative changes, both interactions produce a 
qualitatively similar solution set. Dressing the quark- 
gluon vertex, however, produces qualitative changes as 
well. This is consistent with recent studies that have 
highlighted the impact of vertex dressing on hadron phe- 
nomena [24], HE Hj] an d can be elucidated in the present 
context by charting the solution domains. 

In Fig.[TT]we display the gap equation solution domains 
for the rainbow truncation, whereas those for the 1BC 
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FIG. 11. Chart of gap equation solution domains in the 
{m^,X = D/t<j 2 }-plane obtained with the models specified 
by Eqs. (fit))), ([TTfl . l[T2a)l : top panel - Eq. ifTTajl : and bottom 
panel - Eq. (Illb|) . In both panels the annotations within the 
bounded regions indicate which solutions are found. 



FIG. 12. Chart of gap equation solution domains in the 
{m^,X — D/u; 2 }-plane obtained with the models specified 
by Eqs. ((ID ), l[TT|l . (|12b(l : top panel - Eq. (flla[) : and bottom 
panel - Eq. (jllb|l . In both panels the annotations within the 
bounded regions indicate which solutions are found. 



truncation arc depicted in Fig.[l2j We remark that these 
figures were computed with u) = 0.4 GeV for Eq. (|lla[) 
and uj = 0.5 GeV for Eq. (|llb[) . However, we did vary 
these parameters within the domains described in con- 
nection with Eqs. (fTTj) and found that there is little vari- 
ation. 

It is clear from the figures that little of interest is possi- 
ble until the interaction strength is sufficient to support 
nonperturbative solutions to the gap equation. There- 
after, however, the rainbow truncation can produce a 
novel Wigner-like solution, W2, whose momentum de- 
pendence is typified by Fig. 01 but this is particular to 
that truncation. The studies with a modestly dressed 
vertex show a simpler, regular pattern. We have already 
noted that the gap equation's solution set is quite compli- 
cated at very large interaction strengths. However, such 
strengths far exceed the upper bound on values which are 
capable of describing hadron obscrvables (see footnote [3]) 
and hence we do not describe them herein. 

Let us imagine for the moment that QCD's gap equa- 
tion possesses a kernel whose solution set is one of the 



complicated domains. It may be argued that the different 
solutions within a domain represent competing phases; 
and should they exist simultaneously, then the compu- 
tation of hadron properties would become a complicated 
affair. Moreover, their existence would likely be reflected 
in the properties of hadrons. Since the vast body of 
DSE-based hadron phenomenology does not show any 
sign that this is the case there must be an egress. 

C. Phase Stability 

Egress lies in the direction of phase stability. One must 
consider which of the solutions within a given domain is 
stable against fluctuations. Figures [T3] and [U] contain the 
information necessary to address this question through 
the stability criterion introduced as the generalisation of 
Eqs. ©. 

A careful examination of Fig. [T3] reveals that solutions 
of the inhomogeneous Bethe-Salpetcr equations do not 
exhibit bound-state poles until I > If; i.e., until the in- 
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FIG. 13. X = £> /^-dependence of m 2 an m 2 obtained via 
the Bethe-Salpeter equations in the chiral limit. The pan- 
els depict results obtained with different kernels; namely, 
Eg. (fT0|) and: upper-left - Eqs. (|lla[) . (|12a[) ; upper-right - 
Eqs. (|llb|l , (fT2a|) ; low er-left - Eqs. (fTTa) , ()12b[) ; and lower 
right Eqs. (|llb|) . ()12b|l . All panels: filled squares - m% — rri l a 
along the Wi solution trajectory; open circles - m 2 = m 2 
along Wi\ open up-triangles - m 2 along N^~; filled circles - 
along A"^; open diamonds - m 2 along N^; and down- 



triangles 



along N% 



teraction strength exceeds the critical value for DCSB. 
(This is another example of the causal connection be- 
tween confinement and DCSB in DSE models of QCD 
- see, e.g., Sec. 2 in Ref. Q.) Amidst the solutions dis- 
played beyond X — II, only the regular Nambu solution 
of the gap equation is stable: it produces the well-known 
DCSB case of a masslcss pseudoscalar meson accompa- 
nied by a massive scalar (point "I" in Fig.fJJ). In the chiral 
limit the negative Nambu solution, its partner, produces 
the same results and is equally stable. By the same token, 
the partners to the other displayed solutions are unstable. 
Plainly, apart from the peculiarity of Wi, these results 
are qualitatively the same in all models considered. 

Let us turn now to Fig.fTJ] In all panels the interac- 
tion strength is large enough to induce DCSB at rh = 0; 
and it is abundantly clear that rn^ and rv? a are only 
both positive semi-definite along the trajectory of the 
regular Nambu solution, N± . This is true on an un- 
bounded domain of rh > 0. Unsurprisingly, given the 
qualitative connection between our stability criterion and 
a "Mexican hat" potential, the negative Nambu solution, 
N± , is a saddle-point trajectory: m 2 a > but n\\ < 
(point "S" in Fig.QJ. This nature persists until the 
current-quark mass exceeds a critical value, whose mag- 
nitude is model-dependent but may be characterised as 
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FIG. 14. Current-quark-mass-dependence of m 2 (upper 
grouping) and m 2 , (lower grouping) obtained with different 
kernels; namely, Eq. (|10[) and, within each grouping: upper- 
left - Eos. (fTTa)) , (TrjaJ, X = 5.5; upper-right - Eqs. ([Tib)! , 
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squares 
monds 
stars 



- m 2 j 

m 2 along N 2 ; 



m 2 



along N^; up-triangles 
right-triangles 



along N x ; dia- 
r along JVa"; and 



m, 



along the sole Wigner trajectory. 



m 



i ~ 0.06 ±0.03 GeV. This was first observed in n%m 



Inspection of the lower row in each grouping reveals the 
role of the trajectory as a surrogate Wigner solution 
within a connected, bounded domain of current-quark 
mass, just as was discussed in connection with Fig. [6] It 
is evident from Fig.[14]that all other solutions correspond 
to unstable trajectories. 
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dependence of these and related quantities through the 
temperature induced deconfmement and chiral symmetry 
restoration phase transitions. Notably, here as there, 
vanishes because at 1 < If the pseudoscalar correlation 
involving deconfined quarks possesses neither pseudovec- 
tor nor pseudotensor components; i.e., F = = G = H 
in Eq. (|5c). 

With the introduction of a small current-quark mass 
the transitions become a crossover, as evident in the lower 
grouping of panels in Fig.fJjjJ Nevertheless, deconfme- 
ment is still evident through the sharp rise in the trajec- 
tory associated with below If, a domain whereupon 
the mass-scale determining the magnitude of f^, M(0) 
is seen to switch to the explicit chiral symmetry break- 
ing current-quark mass. These results are qualitatively 
identical to those in Fig. 6 of Ref. (58j . 

The gap equation has long been used as a tool for 
identifying field configurations that may optimally be 
employed in constructing a mean-field approximation, 
or improvements thereof, to a theory's generating func- 
tional. Indeed, truncated gap equations are critical in 
the construction of effective actions for composite oper- 
ators [5^] and therefrom developing models for DCSB 
in hadron physics. In this connection, the gap equa- 
tion solutions have often been interpreted as candidate 
vacua, some of them energetically equivalent but dis- 
tinct, related via global symmetry transformations, in 
the sense first described in connection with the Nambu- 
Jona-Lasinio model (30L |60|. Some external agent, such 
as current-quark mass, then tips the balance in favour 
of one solution, which therefore provides the configura- 
tion around which a model Lagrangian is constructed to 
describe field fluctuations; e.g., Refs. (l9l. l6ll - [63| . Such 
models typically arrive at a potential which expresses 
features that are synonymous with those of the "Mex- 
ican hat". These observations provide another context 
for our results; viz., the models possess far more candi- 
date vacua than practitioners had usually imagined, with 
an hierarchical structure such that, within levels, map- 
pings exist between those solutions related by a symme- 
try transformation. Notwithstanding this, the standard 
Nambu solution of the gap equation is the only one that 
is stable in the presence of a nonzero current-quark mass. 



REMARKS AND SUMMARY 



In Fig. [15] we depict the evolution with interaction 
strength of a range of quantities which typify DCSB in 
hadron physics. Proceeding along the N + trajectory 
from large to small values of X = D/lo 2 in the upper 
grouping of panels, one sees behaviour typical of a Nambu 
to Wigner phase transition. In the chiral limit, all order 
parameters for DCSB vanish at the critical interaction 
strength; and the sharp jump in m w shows coincident 
deconfmement. These panels arc qualitatively identi- 
cal to Fig. 3 in Ref. [58| , which shows the temperature 



We argued, by way of examples, that models of QCD's 
gap equation will typically possess many solutions, a fea- 
ture which owes to the nonlincarity of the equation. Al- 
though simple vertex Ansatze were used, we judge that 
results obtained with the 1BC form exemplify the behav- 
ior one should expect with more realistic models. 

The nature and number of the solutions is readily ex- 
plained and understood. Naturally, in the weak coupling 
limit, only the usual pcrturbative - Wigner type - solu- 
tion is possible. On the other hand, the number of chiral- 
limit solutions evolves with interaction strength, so that 
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at large interaction strengths there are many solutions, 
with distinct pointwise behaviour, that express dynami- 
cal chiral symmetry breaking (DCSB). To be clear: there 
are numerous DCSB solutions in addition to that which 
is usually labelled as the Nambu solution. In response to 
increasing current-quark mass, however, the number of 
solutions decrements uniformly as particular thresholds 
are crossed until, above some value of the mass, only the 
regular Nambu solution remains. 

The gap equation's nonperturbative solutions form a 
hierarchy. In the chiral limit there is a solution within 
each level that preserves chiral symmetry but also a set 
of distinct DCSB solutions that are energetically equiva- 
lent and related by a symmetry transformation. A sym- 
metry transformation does not connect solutions in dif- 
ferent levels, however, and nor are solutions in different 
levels degenerate. 

In the context of composite operator effective actions, 
solutions of the gap equation play the role of candidate 
vacua in the sense that one, from amongst all those avail- 
able, should be chosen as the ground state about which 
dynamical fields may fluctuate. A stability criterion is 
necessary before such a choice can be made. One is 
readily derived from a consideration of the scalar and 
pseudoscalar susceptibilities via their explanation of the 
"Mexican hat" potential and relationship to the fully- 
dressed propagators for composite correlations in the 
scalar and pseudoscalar channels. Fortunately for hadron 
physics phenomcnologics, when applied to the array of 
gap equation solutions, this stability test shows that for 
any nonzero current-quark mass only the regular Nambu 
solution of the gap equation is stable against perturba- 
tions. 
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Appendix A: Homotopy continuation method 

In the context of nonlinear integral equations the ho- 
motopy continuation method (l6| enables one to do more 
than follow a single path to a solution: one can also, e.g., 
switch branches at simple furcation points. The approach 
is therefore more powerful and discriminating than sim- 
ple iteration to a solution. We illustrate aspects of the 
method here using the gap equation as an illustrative 
example. 



To proceed one first converts the integral equation into 
a matrix equation using a discretisation method. The 
Chebyshev expansion scheme described in Ref. [§3] is ef- 
ficient. The gap equation may then be represented as 
follows: 



A; 



B( P f) 



< i < N/2 
N/2 <i<N ' 

where < i, j < N. 



(Al) 



Suppose now that one has a control parameter, A G fl, 
upon which the solution of the gap equation depends. 
Herein A = m, the current-quark mass, or A = D/lo 2 , 
the interaction strength. Given a value of A, the gap 
equation can be understood as the identity 



H(u) = On, 



Xi , un — A . 



(A2) 



where H: R^ -1-1 — > is a smooth mapping on some 
closed domain V € M. N+1 and On is the null-vector in 
M. N . The solutions of Eq. (|A2[) are an inverse image of 
the null- vector. Denoted H~ 1 (0n), in general this inverse 
image describes a collection of smooth curves in R^ -1-1 . 
Importantly, so long as Vw, € V: rk(H'(u)) = N; i.e., the 
derivative has maximal rank throughout T> 1 then each one 
of these curves begins and ends on dT>, the boundary of 
T>, and no two intersect. 

In order to elucidate we will return to interpreting 
A as a control parameter, in which case solutions of 
the gap equation depend parametrically on this vari- 
able: X = A (A). In a typical gap equation study 
one may view the solution process as beginning with 
some small nonnegative value of A, locating the zero; 
then repeating the zero finding steps as A is smoothly 
incremented. With this in mind, suppose that at a 
given value of A = Ai the gap equation has a solution 
X 1 = A(Ai); i.e., F(X 1 ; Ai) = On- Suppose in addition 
that one has already obtained the solution on some do- 
main Ao < A < Ai; i.e., one knows A (A) on this domain, 
and 



Urn det K ' ' ^0, 

A-s-Ai aX 



(A3) 



then A 1 = A(Ai) is readily obtained via straightforward 
iteration from A(Aj~); viz., the solution at some nearby 
A^ < Ai. 

On the other hand, suppose A(Ai) is a solution but 



det 



dF{X-\ x ) 



dX 



0. 



(A4) 



x=x^ 



At such a point A 1 £ R N , either rk(H'(u)) = N or 
rk(H'(u)) ^ N. Consider the first possibility, which cor- 
responds to the curves iJ _1 (0Ar) being smooth. In this 
case 



lim 

A->-Ai 



dX 



dX 



oo 



(A5) 



and X locates a singular point of one of the curves gen- 
erated by H" 1 (On). 
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There are still two possibilities: in the neighbourhood 
of Ai the surface X(X) may either be characterised as 
possessing the form of a straightened S-bend or exhibiting 
a turning point; i.e., bending back upon itself. In the first 
case it might be difficult to obtain the solution at Ai by 
iteration but this straightforward procedure will converge 
at A^ = A + e, where e is a small positive number that 
may be determined empirically. The solution at Ai is 
then bracketed and may be found by interpolation; and 
one can continue the iterative procedure on A > A x . 

The situation is different at a turning point, which, 
in the context of our study, locates the critical current- 
quark mass for the transition between a Nambu solution 
and a Wigner solution. Iteration fails at a turning point. 
In this case one may proceed as follows. Suppose one has 
a solution at A^ : 

X(\i) = {x ,xx,...,xn-i}, / A6 % 
F({x ,xi, . . . ,xn-i}; AT/) = On- 

Now shift Xi —> Xi = Xi + <5, with 5 some small num- 
ber and, typically, i = N/2; hold Xi fixed; and solve by 
iteration the problem 

F({\, x , x n }; Xi) = N . (A7) 

This represents an interchange of roles between the con- 
trol parameter and one element of the solution vector. 
That which was previously the control parameter now 
becomes part of a modified solution vector that is sought 
by iteration. This simple method enables one to join 
and follow the trajectory of the second solution, which 



exists simultaneously on A < Ai with that already ob- 
tained. Once one is sufficiently far removed from Ai on 
this new trajectory, straightforward iteration can again 
be employed. 

Return now to Eq. (|A4[) and consider the remaining 
possibility; in particular, rk(H'(u)) = N — 1. In princi- 
ple this could correspond to one of the curves -ff -1 (0./v) 
terminating within T>. For the gap equation, however, 
this is impossible because it would indicate that there is 
some domain of parameter space in which the gap equa- 
tion does not have a solution. Hence for us this case rep- 
resents a point at which two inverse images of the null 
vector intersect. To be more explicit, we encounter this 
situation when incrementing A = D/ui 2 in the chiral limit 
to arrive at a trifurcation point, whereat a Wigner solu- 
tion connects with a Nambu solution and its reflection. 
At such a location both iteration and the role change 
method of Eq. (|A7[) fail in the sense that neither enables 
the subsequent trajectory of all solutions to be followed. 

To circumvent this difficulty we exploit the current- 
quark mass; i.e., we solve H(u) = A4, where Ai is a 
column vector whose first N/2 elements are zero and the 
next N/2 are m. With careful use of the source term 
provided by the current-quark mass, one eliminates the 
trifurcation point, so that all three solution trajectories 
H~ 1 (Ai) become distinct but remain close. A combina- 
tion of iteration and the role change method can subse- 
quently be used to find and track these solutions. 

We note that in all cases when tracking a solution we 
ensure that Eq. (|A3[) is satisfied at each point X(\) so we 
can be certain that no solution is missed. 
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